Используя явную схему крест и неявную схему, решить начально-краевую задачу для дифференциального уравнения гиперболического типа. Аппроксимацию второго начального условия произвести с первым и со вторым порядком. Осуществить реализацию трех вариантов аппроксимации граничных условий, содержащих производные: двухточечная аппроксимация с первым порядком, трехточечная аппроксимация со вторым порядком, двухточечная аппроксимация со вторым порядком. В различные моменты времени вычислить погрешность численного решения путем сравнения результатов с приведенным в задании аналитическим решением $U(x,t)$. Исследовать зависимость погрешности от сеточных параметров $\tau, h$.

Вариант 7

Начально-краевая задача для дифференциального уравнения гиперболического типа: $$\begin{aligned} & \frac{\partial^2 u}{\partial t^2} + 2 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + 2\frac{\partial u}{\partial x} - 3u\\ &u_x(0,t) = \exp(-t)\cos(2t)\\ &u_x(\frac{\pi}{2},t) = 0\\ &u(x,0) = \exp(-x)\cos x\\ &u_t(x,0) = -\exp(-x)\cos x \end{aligned}$$

Аналитическое решение: $U(x,t) = \exp(-t-x)\cos x\cos(2t)$

## Вариант 10 Начально-краевая задача для дифференциального уравнения гиперболического типа: $$\begin{aligned} & \frac{\partial^2 u}{\partial t^2} + 3 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial u}{\partial x} - u -\cos x\exp(-t)\\ &u_x(0,t) = \exp(-t)\\ &u_x(\pi,t) = -\exp(-t)\\ &u(x,0) = \sin x\\ &u_t(x,0) = -\sin x \end{aligned}$$ Аналитическое решение: $U(x,t) = \exp(-t)\sin x$

Подготовительная часть

К оглавлению

In [1]:
import math
import numpy as np
import plotly
import plotly.graph_objs as go
import matplotlib.pyplot as plt
import time
co = 0
In [2]:
def analytic(x,t):
    return math.exp(-t-x)*math.cos(x)*math.cos(2*t)
def f_1(x,t):
    return math.exp(-t)*math.cos(2*t)
def f_2(x,t):
    return 0
def psi(x,t = None):
    return math.exp(-x)*math.cos(x)
def ksi(x,t = None):
    return -psi(x,t)
        
def MSE(x,t,u):
    error = 0
    for x_ in range(0,len(x)):
        for t_ in range(0,len(t)):
            error += (u[t_][x_] - analytic(x[x_],t[t_]))**2
    return (error/(len(x)*len(t)))/(len(x)*len(t))
    

def normm(x,t,u):
    errors = []
    for i in range(len(u)):
        maximum = 0
        for j in range(len(x)):
            tmp = abs(u[i][j] - analytic(x[j], t[i]))
            if maximum < tmp:
                maximum = tmp
        errors.append(maximum)
    return errors


def save_plot(fig):
  global co
  fig.write_html("plots2/plot" + str(co)+".html")
  co += 1
  print("plot " + str(co) + " saved")


def plot_approximations(u_approx, x):
    approximations = ['2d1', '2d2', '3d2']
    for j in range(len(u_approx)):
        fig = go.Figure()
        fig.update_layout(title = approximations[j][0] + " points " +\
                          approximations[j][-1] + \
                          ("-st order" if approximations[j][-1] == '1' else "-nd order")  +\
                          " approximation",
                         xaxis_title = 'x',
                         yaxis_title = 'U(x,t)', 
                         template = 'plotly_dark')
        for i in range(len(u_approx[j])):
            fig.add_trace(go.Scatter(x = x, y = u_approx[j][i], name = 't = ' + str(i)))
        fig.show()    
        #save_plot(fig)
        
def time_error(x,u,T):
    error = []
    for t in range(len(T)):
        error_tmp = 0
        for i in range(len(x)):
            error_tmp += (u[t][i] - analytic(x[i],T[t]))**2
        #error.append(np.sqrt(error_tmp))
        error.append(error_tmp / len(T))
    return error

def different_time(method,K,N,Tend,approx_f = 'so', approx_type = None,len = 50):
    u = []
    time_step_local = []
    space_step = math.pi/(N-1)
    for K_ in range(K,K+len):
        time_step_local.append(Tend/(K_-1))
        u.append(method(N, K_, approx_f = approx_f)[0])
    return u, time_step_local

def different_h(method, K, N,T, approx_f = 'so', approx_type = None, len = 30, a = 1, b = 2, c = -3):
    global co
    u = []
    X = []
    space_steps_local = []
    for N_ in range(N, 100):
        sigma_tmp = T/(K-1)/(math.pi/(N_ - 1))**2
        if sigma_tmp >= 1:
            print("\033[33m {}" .format('Attention!!! Stable condition:'))
            print('sigma > 1: ', sigma_tmp)
        space_steps_local.append(math.pi/(N_-1))
        X.append([0 + math.pi/2/(N-1)*i for i in range(N)])
        u.append(method(N_, K = K, approx_f = approx_f, approx_type = approx_type, a_ = a, b_ = b, c_ = c)[0])
    return u, space_steps_local, X, 

def space_error(x, u, T):
    error = []
    for x_ in range(len(x)):
        error_tmp = 0
        for t in range(len(T)):
            error_tmp += (u[t][x_] - analytic(x[x_], T[t]))**2
        error.append(error_tmp/len(x))
    return error


def different_space_errors(method, N, K, T, x0 = 0, xl = math.pi, approx_f = 'so',length = 10):
    errors = []
    h = (x0-xl)/(N-1)
    tau = T/(K-1)
    for N_ in range(N,N+length):
        space_step_local = (xl-x0)/(N_-1)
        time_step_local = T/(10*N_)
        x_local = [0 + space_step_local*i for i in range(N_)]
        t_local = [0 + time_step_local * i for i in range(10*N_)]
        #для условия стабильности схемы K = 10*N_
        error = space_error(x_local,method(N = N_, K = 10*N_, T = T, approx_f = approx_f)[0],t_local)
        errors.append(error)
    return errors

def different_time_errors(u, K, N, T, x, t, a = 1, b = 2, c = -3, len = 25):
    time_errors = []
    time_steps_local = []
    for u_, K_ in zip(u, range(K, K+len)):
        time_steps_local.append(T/(K_-1))
        t_local = np.arange(0, T + time_steps_local[-1] - 0.01, time_steps_local[-1])
        time_errors.append(time_error(x,u_,t_local))
    return time_errors

def construct_3d_plot(method, method_name):
    approx_e, X_e, T_e = method(approx_type = '2d2')
    z_plot = []
    for j in range(0, len(X_e), 1):
        tmp = []
        for i in range(0, len(X_e[j]), 1):
            tmp.append(analytic(X_e[j][i], T_e[j][i]))
        z_plot.append(tmp)
    line_marker = dict(color='#0066FF', width=2)
    lines = []
    lines.append(go.Scatter3d(x=X_e[0], y=T_e[0], z=approx_e[0], mode='lines', line=line_marker, name = method_name + ' method', legendgroup=1))
    for i, j, k in zip(X_e[1:], T_e[1:], approx_e[1:]):
        lines.append(go.Scatter3d(x=i, y=j, z=k, mode='lines', line=line_marker, showlegend=False))
    line_marker = dict(color='#ff0000', width=2)
    #my_layout = go.Layout({"title": "Views by publisher", scene = go.Scene(),
    #                       "showlegend": False})
    fig_new = go.Figure(data = lines)
    fig_new.update_layout(title = 'Comparison 3d',template="plotly_dark")
    fig_new.update_scenes(xaxis_title=dict(text = 'x'), yaxis_title = dict(text = 't'), zaxis_title = dict(text = "U(x,t)"))
    fig_new.add_trace(go.Scatter3d(x=X_e[0], y=T_e[0], z=z_plot[0], mode='lines', line=line_marker, name = 'Analytic', legendgroup = 2))
    for i, j, k in zip(X_e[1:], T_e[1:], z_plot[1:]):
        fig_new.add_trace(go.Scatter3d(x=i, y=j, z=k, mode='lines', line=line_marker, showlegend=False))
    fig_new.show()
    #save_plot(fig_new)

    print("Total error = ", norm_error(X_e, T_e, approx_e))
In [3]:
def plot_error(errors, type):
    e_tmp = []
    e_x = []
    for e_ in errors:
        e_x.append(len(e_))
        e_tmp.append(sum(e_)/len(e_))
    fig = go.Figure()
    fig.add_trace(go.Scatter(x = e_x, y = e_tmp))
    fig.update_layout(title = 'Errors with different ' + str(type),
                     xaxis_title = "splits", 
                     yaxis_title = "Error", 
                     template = 'plotly_dark')
    #save_plot(fig)
    fig.show()
    
In [4]:
def fill(space_step, N, K, a, tau, approx_f):
    U = np.zeros([K,N])
    tmp = 0
    for j in range(N):
        U[0][j] = psi(tmp)
        if approx_f == 'fo':
            U[1][j] = psi(tmp) + ksi(tmp)*tau
        if approx_f == 'so':
            U[1][j] = psi(tmp) + ksi(tmp)*tau + a**2*(-math.cos(tmp))*tau
        tmp += space_step
    return U
In [5]:
def plot_analytic(N, T, K, fig = None, x0 = 0, xl = math.pi/2, diff_time = 0):
    u_a = np.zeros((K,N))
    space_step = (xl-x0)/(N-1)
    time_step = T/(K-1)
    x = [0+space_step * i for i in range(N)]
    for k in range(K):
        for j in range(N):
            u_a[k][j] = analytic(j * space_step, k * time_step)
    if fig == None:
        fig = go.Figure()
    fig.add_trace(go.Scatter(x = x, y = u_a[-1], name = 't = ' + str(round(time_step,2))))
    fig.update_layout(title="Analytical Solution",
                      xaxis_title="x",
                      yaxis_title="$U(x,t)$",
                      template = 'plotly_dark')
    fig.show()
    for i in range(1, len(u_a)):
      fig.add_trace(go.Scatter(x = x, y = u_a[i], name = 't = ' + str(round(time_step*(i+1),2))))
    #save_plot(fig)
    fig.show()

Аналитическое решение

К оглавлению

In [6]:
N = 20
T = 10
K = 200
In [7]:
a = 1
b = 2
c = -3
x = [0 + math.pi/2/(N-1)*i for i in range(N)]
t = [0 + T/(K-1)*i for i in range(K)]
In [8]:
plot_analytic(N,T,K)

Явная схема крест

К оглавлению

Теория

В варианте поставлена вторая начально-краевая задача для волнового уравнения. На введенной в случае уравнений параболлической природы пространственно-временной сетке аппроксимируем дифференциальное уравнение следующей схемой:

$$\frac{u_j^{k+1} - 2u_j^k + u_j^{k-1}}{\tau^2} = a^2\frac{u_{j+1}^k -2u_j^k + u_{j-1}^k}{h^2} + O(\tau^2+h^2), \hspace{0.5cm} j = \overline{1,N-1}$$

С шаблоном вида

diagram-20211024 (2).png

Явная конечно-разностная схема условно устойчива с условием $\sigma = \frac{a^2\tau^2}{h^2} < 1$. Это накладывает ограничения на сеточные характеристики.

Для определения $u_j^1$ можно воспользоваться аппроксимацией второго начального условия $$\frac{u_j^1 - u_j^0}{\tau} = \psi_2(x_j) = \xi$$

Для повышения порядка аппроксимации воспользуемся разложением в ряд Тейлора на точном решении в окрестности $t = 0$ и подставим вместо второй производной исходным дифференциальным уравнением. В итоге, имеем

$$u_j^1 = \psi_1(x_j) +psi_2(x_j)\tau + a^2\psi''_1(x_j)\frac{\tau}{2}$$

Реализация

In [9]:
def explicit(N = 20, K = 200, T = 10, approx_f = 'so', a_ = 1, b_ = 2, c_ = -3, approx_type = None):
    space_step = (math.pi/2)/(N-1)
    time_step = T/K
    x = [0 + space_step*i for i in range(N)]
    t = [0 + time_step*i for i in range(K)]
    u = fill(space_step, N, K, a_, time_step, approx_f)
    sigma = a_**2 * time_step**2 / space_step**2
    d = 2
    T_ = []
    X = []
    X = [x, x]
    T_.append([0.0 for _ in x])
    T_.append([time_step for _ in x])
    for k in range(1, K-1):
        for j in range(1, N - 1):
            u[k + 1][j] = \
                u[k][j + 1] *(sigma + time_step**2*b_/(2*space_step))/(d*time_step/2+1) +\
                u[k][j] * (2 - 2*sigma + c_*time_step**2) /(d*time_step/2+1) + \
                u[k][j - 1] *(sigma - time_step**2*b_/(2*space_step))/(d*time_step/2+1) + \
                (u[k - 1][j] * (-1)*(1 - d*time_step/2))/(d*time_step/2+1)
        u[k+1][0] = f_1(None, time_step*k)
        u[k+1][-1] = f_2(None, time_step*k)
        X.append(x)
        T_.append([t[k] for _ in x])
    return u, X, T_

Исследования

In [10]:
u_fo = explicit(N, K, approx_f = 'fo')[0]
u_so = explicit(N, K, approx_f = 'so')[0]
u = [u_fo, u_so]
In [11]:
def plotTime(x, u, t, approx_f):
    fig = go.Figure()
    fig.update_layout(title = "Solution with different time by" + approx_f,
                      template = 'plotly_dark')
    for j in range(2,len(u)):
        fig.add_trace(go.Scatter(x = x, y = u[j], name = "t = " + str(round(t[j], 2))))
    fig.show()
    #save_plot(fig)
In [12]:
plotTime(x,u_fo,t, ' first order initial values approximation')
print("MSE = ", MSE(x, t, u_fo))
plotTime(x,u_so,t, ' second order initial values approximation')
print("MSE = ", MSE(x, t, u_so))
MSE =  1.7276710090483947e-08
MSE =  2.6874675027700173e-07
In [13]:
er = normm(x,t,u_fo)
fig_check_tmp = go.Figure()
fig_check_tmp.add_trace(go.Scatter(x = t, y = er))
fig_check_tmp.update_layout(title = "Error",
                            xaxis_title = "t",
                            yaxis_title = r'$\max |u - analytic|$', 
                            template='plotly_dark')
fig_check_tmp.show()
#save_plot(fig_check_tmp)
In [14]:
def norm_error(x,t,u):
    ans = 0.0
    for i in range(len(u)):
        for j in range(len(u[i])):
            ans += (u[i][j] - analytic(x[i][j], t[i][j]))**2
    return ans/(len(u[0])*len(u))
In [15]:
for approx_f in ['so', 'fo']:
    tmp_u, tmp_t = different_time(explicit, K, N, T, approx_f = approx_f)
    plot_error(different_time_errors(tmp_u, K, N, T, x, t, len = 40), 'time by ' + ("second order" if approx_f[0] == 'f' else "first order") +\
                          " approximation initial values")
    #different_time(explicit,K,N,T,approx_f = 'fo',approx_type = None)

При увеличении количества разбиений по времени (уменьшении шага по времени) происходит следующий раскардаш:

  • В случае аппроксимации начальных условий первого и второго рода виден тренд на уменьшение ошибки с уменьшением шага по времени
In [16]:
errors = []
for approx_f in ['so', 'fo']:
    errors = []
    tmp_u, tmp_h, tmp_x = different_h(explicit, K, N, T, approx_f = approx_f)
    for u_,x_ in zip(tmp_u, tmp_x):
        tmp = space_error(x_,u_,t)
        errors.append(sum(tmp)/len(tmp))
    #different_time(explicit,K,N,T,approx_f = 'fo',approx_type = None)
fig_check = go.Figure()
fig_check.add_trace(go.Scatter(x = tmp_h, y = errors))
fig_check.update_layout(title = 'Errors with different space steps',
                        xaxis_title = 'Space step',
                        yaxis_title = 'Error',
                        template = 'plotly_dark')
fig_check.show()
#save_plot(fig_check)
 Attention!!! Stable condition:
sigma > 1:  1.8380375525067307
 Attention!!! Stable condition:
sigma > 1:  2.036606706378649
 Attention!!! Stable condition:
sigma > 1:  2.2453588937824605
 Attention!!! Stable condition:
sigma > 1:  2.464294114718165
 Attention!!! Stable condition:
sigma > 1:  2.693412369185763
 Attention!!! Stable condition:
sigma > 1:  2.932713657185255
 Attention!!! Stable condition:
sigma > 1:  3.182197978716639
 Attention!!! Stable condition:
sigma > 1:  3.441865333779917
 Attention!!! Stable condition:
sigma > 1:  3.7117157223750876
 Attention!!! Stable condition:
sigma > 1:  3.991749144502152
 Attention!!! Stable condition:
sigma > 1:  4.281965600161109
 Attention!!! Stable condition:
sigma > 1:  4.58236508935196
 Attention!!! Stable condition:
sigma > 1:  4.892947612074704
 Attention!!! Stable condition:
sigma > 1:  5.213713168329341
 Attention!!! Stable condition:
sigma > 1:  5.544661758115872
 Attention!!! Stable condition:
sigma > 1:  5.885793381434295
 Attention!!! Stable condition:
sigma > 1:  6.237108038284612
 Attention!!! Stable condition:
sigma > 1:  6.598605728666823
 Attention!!! Stable condition:
sigma > 1:  6.970286452580926
 Attention!!! Stable condition:
sigma > 1:  7.352150210026923
 Attention!!! Stable condition:
sigma > 1:  7.744197001004813
 Attention!!! Stable condition:
sigma > 1:  8.146426825514595
 Attention!!! Stable condition:
sigma > 1:  8.558839683556272
 Attention!!! Stable condition:
sigma > 1:  8.981435575129842
 Attention!!! Stable condition:
sigma > 1:  9.414214500235303
 Attention!!! Stable condition:
sigma > 1:  9.85717645887266
 Attention!!! Stable condition:
sigma > 1:  10.31032145104191
 Attention!!! Stable condition:
sigma > 1:  10.773649476743053
 Attention!!! Stable condition:
sigma > 1:  11.247160535976088
 Attention!!! Stable condition:
sigma > 1:  11.73085462874102
 Attention!!! Stable condition:
sigma > 1:  12.224731755037842
 Attention!!! Stable condition:
sigma > 1:  12.728791914866555
 Attention!!! Stable condition:
sigma > 1:  13.243035108227165
 Attention!!! Stable condition:
sigma > 1:  13.767461335119668
 Attention!!! Stable condition:
sigma > 1:  14.30207059554406
 Attention!!! Stable condition:
sigma > 1:  14.84686288950035
 Attention!!! Stable condition:
sigma > 1:  15.401838216988534
 Attention!!! Stable condition:
sigma > 1:  15.966996578008608
 Attention!!! Stable condition:
sigma > 1:  16.542337972560578
 Attention!!! Stable condition:
sigma > 1:  17.127862400644435
 Attention!!! Stable condition:
sigma > 1:  17.72356986226019
 Attention!!! Stable condition:
sigma > 1:  18.32946035740784
 Attention!!! Stable condition:
sigma > 1:  18.945533886087386
 Attention!!! Stable condition:
sigma > 1:  19.571790448298817
 Attention!!! Stable condition:
sigma > 1:  20.208230044042143
 Attention!!! Stable condition:
sigma > 1:  20.854852673317364
 Attention!!! Stable condition:
sigma > 1:  21.51165833612448
 Attention!!! Stable condition:
sigma > 1:  22.178647032463488
 Attention!!! Stable condition:
sigma > 1:  22.855818762334387
 Attention!!! Stable condition:
sigma > 1:  23.54317352573718
 Attention!!! Stable condition:
sigma > 1:  24.24071132267187
 Attention!!! Stable condition:
sigma > 1:  24.948432153138448
 Attention!!! Stable condition:
sigma > 1:  25.66633601713692
 Attention!!! Stable condition:
sigma > 1:  26.39442291466729
 Attention!!! Stable condition:
sigma > 1:  27.132692845729547
 Attention!!! Stable condition:
sigma > 1:  27.881145810323705
 Attention!!! Stable condition:
sigma > 1:  28.639781808449754
 Attention!!! Stable condition:
sigma > 1:  29.40860084010769
 Attention!!! Stable condition:
sigma > 1:  30.187602905297528
 Attention!!! Stable condition:
sigma > 1:  30.97678800401925
 Attention!!! Stable condition:
sigma > 1:  31.776156136272867
 Attention!!! Stable condition:
sigma > 1:  32.58570730205838
 Attention!!! Stable condition:
sigma > 1:  33.40544150137579
 Attention!!! Stable condition:
sigma > 1:  34.23535873422509
 Attention!!! Stable condition:
sigma > 1:  35.07545900060627
 Attention!!! Stable condition:
sigma > 1:  35.92574230051937
 Attention!!! Stable condition:
sigma > 1:  36.78620863396435
 Attention!!! Stable condition:
sigma > 1:  37.65685800094121
 Attention!!! Stable condition:
sigma > 1:  38.537690401449986
 Attention!!! Stable condition:
sigma > 1:  39.42870583549064
 Attention!!! Stable condition:
sigma > 1:  40.329904303063195
 Attention!!! Stable condition:
sigma > 1:  41.24128580416764
 Attention!!! Stable condition:
sigma > 1:  42.16285033880398
 Attention!!! Stable condition:
sigma > 1:  43.09459790697221
 Attention!!! Stable condition:
sigma > 1:  44.03652850867233
 Attention!!! Stable condition:
sigma > 1:  44.98864214390435
 Attention!!! Stable condition:
sigma > 1:  45.950938812668255
 Attention!!! Stable condition:
sigma > 1:  46.92341851496408
 Attention!!! Stable condition:
sigma > 1:  47.90608125079176
 Attention!!! Stable condition:
sigma > 1:  48.89892702015137
<ipython-input-2-bd4119510d0c>:94: RuntimeWarning:

overflow encountered in double_scalars

 Attention!!! Stable condition:
sigma > 1:  1.8380375525067307
 Attention!!! Stable condition:
sigma > 1:  2.036606706378649
 Attention!!! Stable condition:
sigma > 1:  2.2453588937824605
 Attention!!! Stable condition:
sigma > 1:  2.464294114718165
 Attention!!! Stable condition:
sigma > 1:  2.693412369185763
 Attention!!! Stable condition:
sigma > 1:  2.932713657185255
 Attention!!! Stable condition:
sigma > 1:  3.182197978716639
 Attention!!! Stable condition:
sigma > 1:  3.441865333779917
 Attention!!! Stable condition:
sigma > 1:  3.7117157223750876
 Attention!!! Stable condition:
sigma > 1:  3.991749144502152
 Attention!!! Stable condition:
sigma > 1:  4.281965600161109
 Attention!!! Stable condition:
sigma > 1:  4.58236508935196
 Attention!!! Stable condition:
sigma > 1:  4.892947612074704
 Attention!!! Stable condition:
sigma > 1:  5.213713168329341
 Attention!!! Stable condition:
sigma > 1:  5.544661758115872
 Attention!!! Stable condition:
sigma > 1:  5.885793381434295
 Attention!!! Stable condition:
sigma > 1:  6.237108038284612
 Attention!!! Stable condition:
sigma > 1:  6.598605728666823
 Attention!!! Stable condition:
sigma > 1:  6.970286452580926
 Attention!!! Stable condition:
sigma > 1:  7.352150210026923
 Attention!!! Stable condition:
sigma > 1:  7.744197001004813
 Attention!!! Stable condition:
sigma > 1:  8.146426825514595
 Attention!!! Stable condition:
sigma > 1:  8.558839683556272
 Attention!!! Stable condition:
sigma > 1:  8.981435575129842
 Attention!!! Stable condition:
sigma > 1:  9.414214500235303
 Attention!!! Stable condition:
sigma > 1:  9.85717645887266
 Attention!!! Stable condition:
sigma > 1:  10.31032145104191
 Attention!!! Stable condition:
sigma > 1:  10.773649476743053
 Attention!!! Stable condition:
sigma > 1:  11.247160535976088
 Attention!!! Stable condition:
sigma > 1:  11.73085462874102
 Attention!!! Stable condition:
sigma > 1:  12.224731755037842
 Attention!!! Stable condition:
sigma > 1:  12.728791914866555
 Attention!!! Stable condition:
sigma > 1:  13.243035108227165
 Attention!!! Stable condition:
sigma > 1:  13.767461335119668
 Attention!!! Stable condition:
sigma > 1:  14.30207059554406
 Attention!!! Stable condition:
sigma > 1:  14.84686288950035
 Attention!!! Stable condition:
sigma > 1:  15.401838216988534
 Attention!!! Stable condition:
sigma > 1:  15.966996578008608
 Attention!!! Stable condition:
sigma > 1:  16.542337972560578
 Attention!!! Stable condition:
sigma > 1:  17.127862400644435
 Attention!!! Stable condition:
sigma > 1:  17.72356986226019
 Attention!!! Stable condition:
sigma > 1:  18.32946035740784
 Attention!!! Stable condition:
sigma > 1:  18.945533886087386
 Attention!!! Stable condition:
sigma > 1:  19.571790448298817
 Attention!!! Stable condition:
sigma > 1:  20.208230044042143
 Attention!!! Stable condition:
sigma > 1:  20.854852673317364
 Attention!!! Stable condition:
sigma > 1:  21.51165833612448
 Attention!!! Stable condition:
sigma > 1:  22.178647032463488
 Attention!!! Stable condition:
sigma > 1:  22.855818762334387
 Attention!!! Stable condition:
sigma > 1:  23.54317352573718
 Attention!!! Stable condition:
sigma > 1:  24.24071132267187
 Attention!!! Stable condition:
sigma > 1:  24.948432153138448
 Attention!!! Stable condition:
sigma > 1:  25.66633601713692
 Attention!!! Stable condition:
sigma > 1:  26.39442291466729
 Attention!!! Stable condition:
sigma > 1:  27.132692845729547
 Attention!!! Stable condition:
sigma > 1:  27.881145810323705
 Attention!!! Stable condition:
sigma > 1:  28.639781808449754
 Attention!!! Stable condition:
sigma > 1:  29.40860084010769
 Attention!!! Stable condition:
sigma > 1:  30.187602905297528
 Attention!!! Stable condition:
sigma > 1:  30.97678800401925
 Attention!!! Stable condition:
sigma > 1:  31.776156136272867
 Attention!!! Stable condition:
sigma > 1:  32.58570730205838
 Attention!!! Stable condition:
sigma > 1:  33.40544150137579
 Attention!!! Stable condition:
sigma > 1:  34.23535873422509
 Attention!!! Stable condition:
sigma > 1:  35.07545900060627
 Attention!!! Stable condition:
sigma > 1:  35.92574230051937
 Attention!!! Stable condition:
sigma > 1:  36.78620863396435
 Attention!!! Stable condition:
sigma > 1:  37.65685800094121
 Attention!!! Stable condition:
sigma > 1:  38.537690401449986
 Attention!!! Stable condition:
sigma > 1:  39.42870583549064
 Attention!!! Stable condition:
sigma > 1:  40.329904303063195
 Attention!!! Stable condition:
sigma > 1:  41.24128580416764
 Attention!!! Stable condition:
sigma > 1:  42.16285033880398
 Attention!!! Stable condition:
sigma > 1:  43.09459790697221
 Attention!!! Stable condition:
sigma > 1:  44.03652850867233
 Attention!!! Stable condition:
sigma > 1:  44.98864214390435
 Attention!!! Stable condition:
sigma > 1:  45.950938812668255
 Attention!!! Stable condition:
sigma > 1:  46.92341851496408
 Attention!!! Stable condition:
sigma > 1:  47.90608125079176
 Attention!!! Stable condition:
sigma > 1:  48.89892702015137
  • В случае невыполнения условия стабильности $\frac{a^2\tau}{h^2}$ решение расходится. В обратном случае -- начинает сходиться.
In [17]:
plot_error(different_space_errors(explicit, N, K, T, approx_f = 'fo'), 'space (first order approximation initial values)')
plot_error(different_space_errors(explicit, N, K, T, approx_f = 'so'), 'space (second order approximation initial values)')
  • При уменьшении шага по пространству (увеличении количества разбиений по пространству) оба типа аппроксимации начальных условий уменьшаются (в случае соблюдения условия стабильности)
In [18]:
construct_3d_plot(explicit, 'explicit')
Total error =  0.0015470497823632438

Теория

Неявная схема записывается в виде:

$$\frac{u_j^{k+1} - 2u_j^k + u_j^{k-1}}{\tau^2} = a^2\frac{u_{j+1}^{k+1} - 2u_j^{k+1} + u_{j-1}^{k+1}}{h^2} + O(\tau + h^2), \hspace{0.5cm} j = \overline{1,N-1}$$

Обладает абсолютной устойчивостью. Сводится к СЛАУ с трехдиагональной матрицей, решаемой методом прогонки. Аппроксимация начальных условий происходит тем же способом, что и в случае явной схемы.

Шаблон схемы

diagram-20211024 (1).png

Реализация

In [19]:
def tma(a,b,c,d,s):
    P = np.zeros(s)
    Q = np.zeros(s)

    P[0] = -c[0] / b[0]
    Q[0] = d[0] / b[0]

    k = s - 1

    for i in range(1, s):
        P[i] = -c[i] / (b[i] + a[i] * P[i - 1])
        Q[i] = (d[i] - a[i] * Q[i - 1]) / (b[i] + a[i] * P[i - 1])
    P[k] = 0
    Q[k] = (d[k] - a[k] * Q[k - 1]) / (b[k] + a[k] * P[k - 1])

    x = np.zeros(s)
    x[k] = Q[k]

    for i in range(s - 2, -1, -1):
        x[i] = P[i] * x[i + 1] + Q[i]
    return x
In [20]:
def implicit(N = 20, K = 200, T = 10,x0 = 0, approx_f = 'so', approx_type = None, a_ = 1, b_ = 2, c_ = -3):
    space_step = (math.pi/2)/(N-1)
    time_step = T/K
    u = fill(space_step, N, K, a_, time_step, approx_f)
    x = [0 + space_step*i for i in range(N)]
    t = [0 + time_step*i for i in range(K)]
    sigma = a_**2 * time_step**2 / space_step**2
    X = []
    T_ = []
    X = [x, x]
    T_.append([0.0 for _ in x])
    T_.append([time_step for _ in x])
    alpha = 0
    beta = 1
    gamma = 0
    delta = 1
    tmp = b_*time_step**2/(2*space_step)
    d_ = 2
    for k in range(1, K - 1):
        a = np.zeros(N)
        b = np.zeros(N)
        c = np.zeros(N)
        d = np.zeros(N)

        for j in range(1, N - 1):
            a[j] = -sigma + tmp
            b[j] = (d_*time_step/2 + 1 + 2*sigma - c_*time_step**2)
            c[j] = -sigma - tmp
            d[j] = u[k - 1][j]*(d_*time_step/2 - 1) + 2*u[k][j]

        # граничные условия не содержат производных
        b[0] = beta - alpha / space_step
        c[0] = alpha / space_step
        d[0] = f_1(None, (k + 1) * time_step)

        a[-1] = - gamma / space_step
        b[-1] = delta + gamma / space_step
        d[-1] = f_2(None, (k + 1) * time_step)

        X.append(x)
        T_.append([t[k] for _ in x])
        Y = tma(a, b, c, d, N)
        u[k + 1] = Y
    return u, X, T_

Исследования

In [21]:
u_fo = implicit(N, K, x0=0, approx_f = 'fo')[0]
u_so = implicit(N, K, x0=0, approx_f = 'so')[0]
u = [u_fo, u_so]
In [22]:
plotTime(x,u_fo,t, ' first order initial values approximation')
print(MSE(x, t, u_fo))
plotTime(x,u_so,t, ' second order initial values approximation')
print(MSE(x, t, u_so))
7.0483381981827445e-09
2.1822005887071565e-07
In [23]:
er = normm(x,t,u_so)
fig_check_tmp = go.Figure()
fig_check_tmp.add_trace(go.Scatter(x = t, y = er))
fig_check_tmp.update_layout(title = "Error",
                            xaxis_title = "t",
                            yaxis_title = r'$\max |u - analytic|$', 
                            template='plotly_dark')
fig_check_tmp.show()
#save_plot(fig_check_tmp)
In [24]:
for approx_f in ['so', 'fo']:
    tmp_u = different_time(implicit, K, N, T, approx_f = approx_f)[0]
    plot_error(different_time_errors(tmp_u, K, N, T, x, t, len = 40), 'time by ' + ("second order" if approx_f[0] == 'f' else "first order") +\
                          " approximation initial values")
    #different_time(explicit,K,N,T,approx_f = 'fo',approx_type = None)
In [25]:
errors = []
for approx_f in ['so', 'fo']:
    errors = []
    tmp_u, tmp_h, tmp_x = different_h(implicit, K, N, T, approx_f = approx_f)
    for u_,x_ in zip(tmp_u, tmp_x):
        tmp = space_error(x_,u_,t)
        errors.append(sum(tmp)/len(tmp))
    #different_time(explicit,K,N,T,approx_f = 'fo',approx_type = None)
fig_check = go.Figure()
fig_check.add_trace(go.Scatter(x = tmp_h, y = errors))
fig_check.update_layout(title = 'Error with different space steps',
                        xaxis_title = 'Space step',
                        yaxis_title = 'Error',
                        template = 'plotly_dark')
fig_check.show()
 Attention!!! Stable condition:
sigma > 1:  1.8380375525067307
 Attention!!! Stable condition:
sigma > 1:  2.036606706378649
 Attention!!! Stable condition:
sigma > 1:  2.2453588937824605
 Attention!!! Stable condition:
sigma > 1:  2.464294114718165
 Attention!!! Stable condition:
sigma > 1:  2.693412369185763
 Attention!!! Stable condition:
sigma > 1:  2.932713657185255
 Attention!!! Stable condition:
sigma > 1:  3.182197978716639
 Attention!!! Stable condition:
sigma > 1:  3.441865333779917
 Attention!!! Stable condition:
sigma > 1:  3.7117157223750876
 Attention!!! Stable condition:
sigma > 1:  3.991749144502152
 Attention!!! Stable condition:
sigma > 1:  4.281965600161109
 Attention!!! Stable condition:
sigma > 1:  4.58236508935196
 Attention!!! Stable condition:
sigma > 1:  4.892947612074704
 Attention!!! Stable condition:
sigma > 1:  5.213713168329341
 Attention!!! Stable condition:
sigma > 1:  5.544661758115872
 Attention!!! Stable condition:
sigma > 1:  5.885793381434295
 Attention!!! Stable condition:
sigma > 1:  6.237108038284612
 Attention!!! Stable condition:
sigma > 1:  6.598605728666823
 Attention!!! Stable condition:
sigma > 1:  6.970286452580926
 Attention!!! Stable condition:
sigma > 1:  7.352150210026923
 Attention!!! Stable condition:
sigma > 1:  7.744197001004813
 Attention!!! Stable condition:
sigma > 1:  8.146426825514595
 Attention!!! Stable condition:
sigma > 1:  8.558839683556272
 Attention!!! Stable condition:
sigma > 1:  8.981435575129842
 Attention!!! Stable condition:
sigma > 1:  9.414214500235303
 Attention!!! Stable condition:
sigma > 1:  9.85717645887266
 Attention!!! Stable condition:
sigma > 1:  10.31032145104191
 Attention!!! Stable condition:
sigma > 1:  10.773649476743053
 Attention!!! Stable condition:
sigma > 1:  11.247160535976088
 Attention!!! Stable condition:
sigma > 1:  11.73085462874102
 Attention!!! Stable condition:
sigma > 1:  12.224731755037842
 Attention!!! Stable condition:
sigma > 1:  12.728791914866555
 Attention!!! Stable condition:
sigma > 1:  13.243035108227165
 Attention!!! Stable condition:
sigma > 1:  13.767461335119668
 Attention!!! Stable condition:
sigma > 1:  14.30207059554406
 Attention!!! Stable condition:
sigma > 1:  14.84686288950035
 Attention!!! Stable condition:
sigma > 1:  15.401838216988534
 Attention!!! Stable condition:
sigma > 1:  15.966996578008608
 Attention!!! Stable condition:
sigma > 1:  16.542337972560578
 Attention!!! Stable condition:
sigma > 1:  17.127862400644435
 Attention!!! Stable condition:
sigma > 1:  17.72356986226019
 Attention!!! Stable condition:
sigma > 1:  18.32946035740784
 Attention!!! Stable condition:
sigma > 1:  18.945533886087386
 Attention!!! Stable condition:
sigma > 1:  19.571790448298817
 Attention!!! Stable condition:
sigma > 1:  20.208230044042143
 Attention!!! Stable condition:
sigma > 1:  20.854852673317364
 Attention!!! Stable condition:
sigma > 1:  21.51165833612448
 Attention!!! Stable condition:
sigma > 1:  22.178647032463488
 Attention!!! Stable condition:
sigma > 1:  22.855818762334387
 Attention!!! Stable condition:
sigma > 1:  23.54317352573718
 Attention!!! Stable condition:
sigma > 1:  24.24071132267187
 Attention!!! Stable condition:
sigma > 1:  24.948432153138448
 Attention!!! Stable condition:
sigma > 1:  25.66633601713692
 Attention!!! Stable condition:
sigma > 1:  26.39442291466729
 Attention!!! Stable condition:
sigma > 1:  27.132692845729547
 Attention!!! Stable condition:
sigma > 1:  27.881145810323705
 Attention!!! Stable condition:
sigma > 1:  28.639781808449754
 Attention!!! Stable condition:
sigma > 1:  29.40860084010769
 Attention!!! Stable condition:
sigma > 1:  30.187602905297528
 Attention!!! Stable condition:
sigma > 1:  30.97678800401925
 Attention!!! Stable condition:
sigma > 1:  31.776156136272867
 Attention!!! Stable condition:
sigma > 1:  32.58570730205838
 Attention!!! Stable condition:
sigma > 1:  33.40544150137579
 Attention!!! Stable condition:
sigma > 1:  34.23535873422509
 Attention!!! Stable condition:
sigma > 1:  35.07545900060627
 Attention!!! Stable condition:
sigma > 1:  35.92574230051937
 Attention!!! Stable condition:
sigma > 1:  36.78620863396435
 Attention!!! Stable condition:
sigma > 1:  37.65685800094121
 Attention!!! Stable condition:
sigma > 1:  38.537690401449986
 Attention!!! Stable condition:
sigma > 1:  39.42870583549064
 Attention!!! Stable condition:
sigma > 1:  40.329904303063195
 Attention!!! Stable condition:
sigma > 1:  41.24128580416764
 Attention!!! Stable condition:
sigma > 1:  42.16285033880398
 Attention!!! Stable condition:
sigma > 1:  43.09459790697221
 Attention!!! Stable condition:
sigma > 1:  44.03652850867233
 Attention!!! Stable condition:
sigma > 1:  44.98864214390435
 Attention!!! Stable condition:
sigma > 1:  45.950938812668255
 Attention!!! Stable condition:
sigma > 1:  46.92341851496408
 Attention!!! Stable condition:
sigma > 1:  47.90608125079176
 Attention!!! Stable condition:
sigma > 1:  48.89892702015137
 Attention!!! Stable condition:
sigma > 1:  1.8380375525067307
 Attention!!! Stable condition:
sigma > 1:  2.036606706378649
 Attention!!! Stable condition:
sigma > 1:  2.2453588937824605
 Attention!!! Stable condition:
sigma > 1:  2.464294114718165
 Attention!!! Stable condition:
sigma > 1:  2.693412369185763
 Attention!!! Stable condition:
sigma > 1:  2.932713657185255
 Attention!!! Stable condition:
sigma > 1:  3.182197978716639
 Attention!!! Stable condition:
sigma > 1:  3.441865333779917
 Attention!!! Stable condition:
sigma > 1:  3.7117157223750876
 Attention!!! Stable condition:
sigma > 1:  3.991749144502152
 Attention!!! Stable condition:
sigma > 1:  4.281965600161109
 Attention!!! Stable condition:
sigma > 1:  4.58236508935196
 Attention!!! Stable condition:
sigma > 1:  4.892947612074704
 Attention!!! Stable condition:
sigma > 1:  5.213713168329341
 Attention!!! Stable condition:
sigma > 1:  5.544661758115872
 Attention!!! Stable condition:
sigma > 1:  5.885793381434295
 Attention!!! Stable condition:
sigma > 1:  6.237108038284612
 Attention!!! Stable condition:
sigma > 1:  6.598605728666823
 Attention!!! Stable condition:
sigma > 1:  6.970286452580926
 Attention!!! Stable condition:
sigma > 1:  7.352150210026923
 Attention!!! Stable condition:
sigma > 1:  7.744197001004813
 Attention!!! Stable condition:
sigma > 1:  8.146426825514595
 Attention!!! Stable condition:
sigma > 1:  8.558839683556272
 Attention!!! Stable condition:
sigma > 1:  8.981435575129842
 Attention!!! Stable condition:
sigma > 1:  9.414214500235303
 Attention!!! Stable condition:
sigma > 1:  9.85717645887266
 Attention!!! Stable condition:
sigma > 1:  10.31032145104191
 Attention!!! Stable condition:
sigma > 1:  10.773649476743053
 Attention!!! Stable condition:
sigma > 1:  11.247160535976088
 Attention!!! Stable condition:
sigma > 1:  11.73085462874102
 Attention!!! Stable condition:
sigma > 1:  12.224731755037842
 Attention!!! Stable condition:
sigma > 1:  12.728791914866555
 Attention!!! Stable condition:
sigma > 1:  13.243035108227165
 Attention!!! Stable condition:
sigma > 1:  13.767461335119668
 Attention!!! Stable condition:
sigma > 1:  14.30207059554406
 Attention!!! Stable condition:
sigma > 1:  14.84686288950035
 Attention!!! Stable condition:
sigma > 1:  15.401838216988534
 Attention!!! Stable condition:
sigma > 1:  15.966996578008608
 Attention!!! Stable condition:
sigma > 1:  16.542337972560578
 Attention!!! Stable condition:
sigma > 1:  17.127862400644435
 Attention!!! Stable condition:
sigma > 1:  17.72356986226019
 Attention!!! Stable condition:
sigma > 1:  18.32946035740784
 Attention!!! Stable condition:
sigma > 1:  18.945533886087386
 Attention!!! Stable condition:
sigma > 1:  19.571790448298817
 Attention!!! Stable condition:
sigma > 1:  20.208230044042143
 Attention!!! Stable condition:
sigma > 1:  20.854852673317364
 Attention!!! Stable condition:
sigma > 1:  21.51165833612448
 Attention!!! Stable condition:
sigma > 1:  22.178647032463488
 Attention!!! Stable condition:
sigma > 1:  22.855818762334387
 Attention!!! Stable condition:
sigma > 1:  23.54317352573718
 Attention!!! Stable condition:
sigma > 1:  24.24071132267187
 Attention!!! Stable condition:
sigma > 1:  24.948432153138448
 Attention!!! Stable condition:
sigma > 1:  25.66633601713692
 Attention!!! Stable condition:
sigma > 1:  26.39442291466729
 Attention!!! Stable condition:
sigma > 1:  27.132692845729547
 Attention!!! Stable condition:
sigma > 1:  27.881145810323705
 Attention!!! Stable condition:
sigma > 1:  28.639781808449754
 Attention!!! Stable condition:
sigma > 1:  29.40860084010769
 Attention!!! Stable condition:
sigma > 1:  30.187602905297528
 Attention!!! Stable condition:
sigma > 1:  30.97678800401925
 Attention!!! Stable condition:
sigma > 1:  31.776156136272867
 Attention!!! Stable condition:
sigma > 1:  32.58570730205838
 Attention!!! Stable condition:
sigma > 1:  33.40544150137579
 Attention!!! Stable condition:
sigma > 1:  34.23535873422509
 Attention!!! Stable condition:
sigma > 1:  35.07545900060627
 Attention!!! Stable condition:
sigma > 1:  35.92574230051937
 Attention!!! Stable condition:
sigma > 1:  36.78620863396435
 Attention!!! Stable condition:
sigma > 1:  37.65685800094121
 Attention!!! Stable condition:
sigma > 1:  38.537690401449986
 Attention!!! Stable condition:
sigma > 1:  39.42870583549064
 Attention!!! Stable condition:
sigma > 1:  40.329904303063195
 Attention!!! Stable condition:
sigma > 1:  41.24128580416764
 Attention!!! Stable condition:
sigma > 1:  42.16285033880398
 Attention!!! Stable condition:
sigma > 1:  43.09459790697221
 Attention!!! Stable condition:
sigma > 1:  44.03652850867233
 Attention!!! Stable condition:
sigma > 1:  44.98864214390435
 Attention!!! Stable condition:
sigma > 1:  45.950938812668255
 Attention!!! Stable condition:
sigma > 1:  46.92341851496408
 Attention!!! Stable condition:
sigma > 1:  47.90608125079176
 Attention!!! Stable condition:
sigma > 1:  48.89892702015137
In [26]:
plot_error(different_space_errors(implicit, N, K, T, approx_f = 'fo'), 'space (first order approximation initial values)')
plot_error(different_space_errors(implicit, N, K, T, approx_f = 'so'), 'space (second order approximation initial values)')
In [27]:
construct_3d_plot(implicit, 'implicit')
Total error =  0.0013915009590317656
In [ ]: